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Abstract 

The Heun functions have wide application in modern physics and are expected 
to succeed the hypergeometrical functions in the physical problems of the 21st 
century. The numerical work with those functions, however, is complicated and 
requires filling the gaps in the theory of the Heun functions and also, creating 
new algorithms able to work with them efficiently. 

We propose a new algorithm for solving a system of two nonlinear transcen- 
dental equations with two complex variables based on the Miiller algorithm. The 
new algorithm is particularly useful in systems featuring the Heun functions and 
for them, the new algorithm gives distinctly better results than Newton's and 
Broydcn's methods. 

As an example for its application in physics, the new algorithm was used 
to find the quasi-normal modes (QNM) of Schwarzschild black hole described 
by the Regge- Wheeler equation. The numerical results obtained by our method 
are compared with the already published QNM frequencies and are found to 
coincide to a great extent with them. Also discussed are the QNM of the Kerr 
black hole, described by the Teukolsky Master equation. 

Keywords: root-finding algorithm, Miiller algorithm, two-dimensional Miiller 
algorithm, Regge- Wheeler equation, QNM, Teukolsky Master equation 
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1. Overview 

The Heun functions appear with increasing frequency in modern physics. 
For example, they arise in the Schrodinger equation with anharmonic potential, 
in water molecule, in the Stark effect, in different quantum phenomena related 
with repulsion and attraction of levels, in the theory of lunar motion, in grav- 
itational physics of scalar, spinor, electromagnetic and gravitational waves in 
Schwarzschild and Kerr metric, in crystalline materials, in three-dimensional 
waves in atmosphere, in Bethe anzatz systems, in Collogero-Moser-Sutherland 
systems, e.t.c, just to mention a few. Because of the wide range of their appli- 
cations ([212])- from quantum mechanics to astrophysics, from lattice systems 
to economics - they can be considered as the 21 st century successors of the 
hypergeometric functions encountered in some simple physical problems of 20 th 
century. 

It is not hard to explain this situation. In natural sciences, in particular, 
in physics we usually study the different phenomena starting from some equi- 
librium state. Then we study small deviations from it in linear approximation, 
and at the end, going far away from the equilibrium we are forced to take into 
account nonlinear phenomena. It is well known that to describe the wave pro- 
cesses (like those in quantum mechanics) , related with some linear phenomenon 
in classical physics (like classical mechanics), we have to use hypergeometric 
functions. Therefore these functions were well studied in 19 th and 20 th centuries 
and today one can find the corresponding codes in all good computer pack- 
ages. According to the theorem by Slavyanov [2J, if we study nonlinear classical 
phenomena, described by elliptic functions, or even by the solutions of any of 
Painleve type equations, the corresponding wave problems can be solved exactly 
in terms of the Heun functions. Since the Painleve equations can be considered 
as Hamilton ones for a very large class of nonlinear classical problems, one can 
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expect a fast increase in the applications of the Heun functions in physics and 
other natural sciences of 21 th century. 

Their mathematical complexity, however, makes working with them a sig- 
nificant challenge both analytically and numerically. The Heun functions are 
unique local Frobenius solutions of a second-order linear ordinary differential 
equation of the Fuchsian type [H [3J HI [5] which in the general case have 4 reg- 
ular singular points. Two or more of those regular singularities can coalesce 
into an irregular singularity leading differential equations of the confluent type 
and their solutions: confluent Heun function, biconfluent Heun function, double 
confluent Heun function and triconfluent Heun function. The Heun functions 
generalize the hypergeometric function (and also include the Lame function, 
Mathieu function and the spheroidal wave functions [H [5]) and some of their 
uses can be found in [2] and also in the more recent [6]. Clearly, the Heun 
functions will be encountered more and more in modern physics, hence, there 
is a need of better understanding of those functions and new, more adequate 
algorithms working with them. 

Despite the growing number of articles which use equations from the Heun 
type and their solutions, the theory of those functions is far from complete. 
There are some analytical works on the Heun functions, but they were largely 
neglected until recently. Some recent progress can be found in [7] , but as a whole 
there are many gaps in our knowledge of those functions. Particularly, the con- 
nection problem for the Heun functions is not solved - one cannot connect two 
local solutions at different singular points using known constant coefficients 0. 
Another example of a serious gap in the general theory of the Heun functions in 
general is the absence of integral representations analogous to the one for hyper- 
geometric functions. According to Whittaker's hypothesis, the Heun's functions 
are the simplest class of special functions for which no representations in form 
of contour integrals of elementary functions exists. In some particular cases 
are known integral representations in form of double integrals of elementary 
functions, but the general case is an open problem. 

Numerically, the only software package currently able to work with the Heun 
functions is MAPLE. Alternative ways for evaluations of those functions do not 
exist (to the best of our knowledge) and there are no known projects aiming 
to change this situation, an admittedly immense task by itself. This means 
that the use of the Heun functions is limited to the routines hidden in the 
kernel of maple, which the user cannot change or improve - a situation that 
makes understanding the numerical problems or avoiding them adequately very 
difficult. On the positive side, those routines were found by the team to work 
well enough in many cases (see, for example, the match between theory and 
numerical results in [8], as well as the other applications of those functions in 
HO] ) . Yet, there are some peculiarities - there are values of the parameters 
where the routines break down leading to infinities or to numerical errors. The 
situation with the derivatives of the Heun functions in maple is even worst 
- for some values they simply do not work, for example outside the domain 
\z\ < I, where their precision is much lower than that of the Heun function itself. 
Also, in some cases there are no convenient power-series representations and 
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then the Heun functions are evaluated in maple using numerical integration. 
Therefore the procedure goes slowly in the complex domain (compared to the 
hypergeometric function) which means that the convergence of the root-finding 
algorithm is essential when one solves equations including Heun's functions. 

Despite all the numerical set-backs, the Heun functions offer many oppor- 
tunities to modern physics. They occur in the problem of quasi-normal modes 
(QNM) of rotating and nonrotating black holes, which is to some extent the 
gravity analogue of the problem of the hydrogen atom. Finding the QNMs is 
critical to understanding observational data from gravitational wave detectors 
and proving or refuting the black holes existence ([TTJ [T2] and also [HI US ITU]). 
In this case, one has to solve a two-dimensional connected spectral problem 
with two complex equations in each of which one encounters the confluent Heun 
functions. The analytical theory of the confluent Heun function is more devel- 
oped than that of the other types of Heun functions, but still many unknowns 
remain. Again, the evaluation of the derivative of the confluent Heun function is 
problematic in maple, which makes Newton's root-finding algorithm ( [131 114) ) 
unusable. Broyden's algorithm ([H]) generally works, but it is slowly convergent 
even close to a root (see [IS]). It is clear that we need a novel algorithm, that 
will offer quicker convergence than Broyden's algorithm, but without relying on 
derivatives. 

To solve this problem, in the case of a system of two equations in two vari- 
ables, our team developed a two-dimensional generalization of the Miiller algo- 
rithm. The one-dimensional Miiller algorithm (|17j) is a quadratic generalization 
of the secant method, that works well in the case of a complex function of one 
variable. It has very good convergence for a large class of functions (~ 1.84) and 
it is very efficient when the starting point (the initial guess) is close to a root 
(for applications see [TS] and [5])- It is also well convergent when working with 
special transcendental functions such as the Heun functions. Most importantly, 
this algorithm does not use derivatives, which is key for our work. 

The two-dimensional Miiller algorithm seems to inherit some of the advan- 
tages of its one-dimensional counterpart like good convergence and usability 
on large class of functions as our tests show. The new algorithm was used to 
solve the QNM problem in the case of a Schwarzschild black hole and it proved 
to work without significant deviations from the results published by Andersson 
(|19j) and Fiziev ([IB]). Also, preliminary results for the QNM of the Kerr black 
hole are discussed and for them we also obtain a very good coincidence with 
published results [2"0] . 

The article is organized as follows: Section 2 reviews the one-dimensional 
Miiller algorithm and it introduces in detail its two-dimensional generalization, 
in Section 3 we give two examples of the application of the new method in 
systems featuring the confluent Heun functions - the QNMs of rotating and 
nonrotating BH, and the results in both cases are analyzed in terms of precision 
and convergence. In Section 4 we summarize our results. 
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2. The Miiller algorithm 



2.1. One- dimensional Mutter's algorithm 

The one-dimensional Miiller algorithm ( [131 117j) is iterative method which 
at each step evaluates the function at three points, builds the parabola crossing 
those points and finds the two points where that parabola crosses the x-axis. 
The next iteration is the the point farthest from the initial point. 

Explicitly, the one-dimensional Miiller algorithm can be defined as the map 
(x : fi(x° , F(x)) — > x p , which obtains the final point x p in P iterations by 
calculating for every three points Xj~2, Xj-i, Xj (with the corresponding values 
of the function f(x) — > fj—i> fj), the next iteration Xj+i as: 

x J+ i = xj - (xj-Xj-i) ^ where 

max{D 1 , D 2 ) 

A = f ] q-q(l + q)f J _ 1 +q 2 f^ 2 ,B = (2q+l)f J -(l + q) 2 f^ 1 +q 2 fj-2, 



C = (l+q)fj, Dus = B± VB 2 -4 AC and 



x j x j-i 
x 3-l ~ x j-2 



The iterations continue until | xp — %p-\ |< 10~ d , where d is the number of 
digits of precision we require. This exit-condition works independently of the 
actual numerical zero in use, which may vary for the confluent Heun function 
and thus it is the most appropriate for our numerical work. 

2.2. Two-dimensional Mutter's algorithm 

The two-dimensional Miiller method comes as a natural extension of the 
one-dimensional Miiller method. 

For two complex- valued functions F\(x, y) and F2(x, y) we want to find such 
pairs of complex numbers (x r , y t ) which are solutions of the system: 

(F 1 (x I ,y I )=0 
\^ 2 (x 7 , yj ) = 

where / = 1, . . . numbers the solution in use. From now on, we will omit 
the index /, considering that we work with one arbitrary particular solution. 
Finding all the solutions of a system is beyond the scope of this article. 

Consider the functions Fi(x,y), F2(x,y) as two-dimensional complex sur- 
faces z = Fi(x,y) and z = F2(x,y) in a three-dimensional space of the complex 
variables {x, y, Normally to solve the system, one expresses the relation 
y(x) from one of the equations, then by substituting it in the other equation, 
one solves it for x and from y(x) one finds y. In the general case, however, this is 
not possible. The idea of our code is to approximately follow that procedure by 
finding an approximate linear relation y(x) between the two variables and then 



1 Equivalently, we can consider four real surfaces in five-dimensional real hyperspace, which 
are defined by four real functions of four real variables 
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A : CiXi + C 2 yi + C 3 = F 2 {xi,yi), i = [n-2, n-l,n] -> Ci,C 2 .C"3 



A*(a:„,Fi(T,t/ = 1/(3-))) -> £„+i 



M2: i-i(y„,F 2 (x = x„ +1 ,y)) -> ij n+1 or Ml: j/ n+1 = y(.r„ +1 ) 



(s^+i.t/n+i) -> /'1 1 , y„ 1 1 ), /-sC 1 " 1 1 , y n ) 1 ) 



Figure 1: A block scheme of the two-dimensional Miiller algorithm. A(Ci, C2, C3) is the plane 
with equation 2 = Cix + C2J/ + C3 that crosses trough the 3 pairs of points (xi,yi) and the 
function F2 evaluated in them. The plane v is defined by the equation 2 = 0. The one- 
dimensional Miiller algorithm, fj,(t ln , F(t)) — > t^ ln , is applied on the function of one variable 
F(t) with starting point t ln and final point t^ In . 



using it to find the root of function of one variable trough the one-dimensional 
Miiller algorithm. 

To find the linear relation y{x), at each iteration we form the plane passing 
trough three points of one of the functions and then the equation of the line of 
intersection between that plane and the plane z = is used as the approximate 
relation y(x). This basically means that the so found y(x) is an approximate 
solution of one of the equations which ideally should be near the real solution in 
the z = plane. Substituting this relation in the other function, we run the one- 
dimensional Miiller algorithm on it to fix the value of one of the variables, say 
x. Using the value of x in the first function, we again run the one-dimensional 
Miiller algorithm on it to fix the value of the other variable - y. Alternatively 
one can substitute the value of x directly in y(x) to obtain y. This ends one 
iteration of the algorithm. The process repeats until one of the exit-conditions 
discussed below has been reached. The block-scheme of the algorithm can be 
seen on Fig. ([!]). 

Explicitly, the code starts by evaluating the two functions Fi^{xi,yi) in 
three starting pairs of points (i = 1,2,3) that ideally should be near one of the 
roots of the system. In our case, those three initial pairs are obtained from 
one starting pair to which we add and subtract certain small complex number. 
This artificial choice is done only in the first iteration (n — 3), afterwards we 
use the output of the last three iterations to form (x n —2, yn—2), (%n—i> J/n-i)j 
{xniVn) an d the respective Fi 2 {x,y). Thus on every iteration after n — 3 the 
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actual complex functions Fi_ 2 (x n ,y n ) are evaluated only once outside of the 
one-dimensional Miiler subroutines. 

Next we construct the plane passing trough those three points for one of the 
functions, say F 2 by solving the linear system: 



From it one obtains the coefficients C\, C 2 , C3 of the plane z = C\x + C' 2 y + C3. 

This plane is intersected with the plane z = and the equation of the line 
between those two planes is the approximate relation y(x) (i.e. y = {—C\x — 
C^)/C 2 ) of the two variables. 

We substitute that relation in the first function Fi(x,y) — > Fi(x,y(x)) and 
we start the one-dimensional Mullcr on that "linearized" function of only one 
variable, x. After some prc-dctcrmincd maximal number of iterations, the exit- 
ing point is chosen for x n+1 ( fi(x n , F^x, y = y(x))) -» x n+1 ) . 

Then, there are two possibilities. 

Algorithm Ml: one could use directly the relation y(x = x n+ i) to find 
y = y n+1 . Or, 

Algorithm M2: One can substitute x — x n+i in the other function F 2 (x, y) — > 
F 2 (x — x n+ i,y) in order to find y n +i using again the one-dimensional Mullcr 
algorithm (n(y n , F 2 (x n+1 , y)) -> y n +i)- 

Our numerical experiments showed that both approaches lead to convergent 
procedure. 

After (x n+ i,y n+ i) are fixed, the two functions F\ t2 (x n +i, y n +i) are evaluated 
and if the new points arc not roots, the iterations continue. 

The exit-strategy in the two-dimensional Muller algorithm is as follows: 

1. To avoid hanging of the algorithm or its deviation from the actual root 
of the system, one fixes the maximal number of iterations for the one- 
dimensional Muller subroutine, P. From our experience, P = 4— 10 gives 
best convergence. 

2. The precision-condition (| Xj — Xj-i |< 10~ d ) remains in force for the 
one-dimensional Muller algorithm. Usually the subroutine exits, because 
of j > P during the first few iterations of the two-dimensional Muller 
algorithm, until it gets closer to the roots and exits with | Xj—Xj-i |< lCr d . 

3. The primary exit-condition of the two-dimensional Mullcr is fulfilled when 
for two consecutive pairs (x n ,y n ): | x n — |< 10~ d , | y n —y n -\ |< 10 _<i . 
In this case, if the values of functions F\{x, y),F 2 (x, y) at those points are 
sufficiently small, the algorithm exits with a root. 

4. To keep the two-dimensional Mullcr algorithm from hanging, a maximal 
number of iterations N should be set. For n > N the algorithm exits 
without fixing a root. 

5. A common problem occurs when one of the functions becomes zero be- 
fore the other function. A possible way out is to substitute one of the 



C\x n - 2 + C 2 y n - 2 + C 3 = F 2 (x 
C\x n -\ + C 2 y n -i + C 3 = F 2 {x n -i,y n -i 
dx n + C 2 y n + C 3 = F 2 (x n ,y n ). 



) 
) 
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so-fixed variables, say x^ ln , in the non-zero function and to to run the 
one-dimensional Muller subroutine using the other variable - /i(y ra , F2(x = 
x i V)) y^ m - The algorithm then exits with a possible root: (x^ m , y^ m ) 

The procedure can be fine-tuned trough change in the starting pair of points, 
the initial deviation or by switching the places of the functions, or even by 
replacing the functions with their independent linear combinations. 

The numerical experiments of the application of this algorithm on systems 
featuring different classes of functions can be found in [16 . The tests showed 
that the algorithm inherits some of the advantages of the one-dimensional Muller 
algorithm, like the quick convergence in proximity of the root and the vast class 
of functions that it can work with. The major disadvantage comes from the 
complicated behavior of the two-dimensional complex surfaces defined by the 
functions ^1.2 (£, y) which requires one to find the best combination of starting 
points and number of iterations in the one-dimensional Muller subroutine so 
that the algorithm converges to the required root (if it is known or suspected). 
Generally, it is hard to tell when one point is "close" to a root. In some cases, 
even if certain starting pair of points is close to a root in terms of some norm, 
using it as a starting point in the algorithm may still lead to convergence to 
another root or simply to require more iterations to reach the desired root than 
if other pair of starting points were used. 

It is important to note that unlike Broyden's algorithm and Newton's algo- 
rithm which do not dependent on the order of the equations in the system, our 
two-dimensional Muller algorithm depends on the order of the equations. The 
numerical experiments show that while for some systems, changing the places 
of the equations has little or no effect on the convergence, in other cases, it 
slows down or completely breaks down the convergence. While such inherent 
asymmetry certainly is a weakness of the algorithm, there are ways around it. 
For example, one may alternate the places of the equations at each iteration 
or use their independent linear combinations (-Fx 2 = c*i,2-Fi + pxpFz). Those 
approaches make the algorithm more robust, but since they may cost speed, we 
prefer to set the order of the equations manually. 

A technical disadvantage is that the whole procedure is more CPU-expensive 
than Newton's method and Broyden's method, since it generally makes more 
evaluations of the functions - each one-dimensional Muller makes at least 1 
iteration on every step of the two-dimensional Muller, thus it makes at least 
4 evaluations of each function. This is because on each iteration of the two- 
dimensional Muller algorithm the functions in use change and thus one cannot 
use previous evaluations to reduce time. Still, in some cases, as demonstrated 
in |16| and also below, the so-constructed algorithm is quicker or comparable to 
Newton's or Broyden's method. 
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3. Some applications of the method QNMs of nonrotating and ro- 
tating black holes 

We will work only with the confluent Heun functions, which are much better 
studied than the other types of Heun functions, due to their numerous physical 
applications. Besides their numerical implementation was used successfully in 
previous works by the authors. For details on the numerical testing, see |16| . 

The quasi-normal modes (QNMs) of a black hole (BH) are the complex 
frequencies (cj) that govern the late-time evolution of the perturbations of the 
BH metric ( [HI [H ED E1I23] ) , which have been extensively studied. 

3.1. First example: Non-rotating black hole 

First, we consider the problem of the gravitational QNMs of a nonrotating 
black hole so that the precision of the new method can be tested on a very well 
studied physical problem. The physical results in this case were published in 
[5], so here we will focus on the numerical details instead. 

To find the QNMs, one uses the exact solutions of the Regge- Wheeler equa- 
tions, describing the linearized perturbations of Schwarzschild metric, in terms 
of confluent Heun functions([18 j ). From [18], when the mass of the BH is set to 
2M = 1, one obtains the following system of the type ([I]): 

Fi = (cos(0) - l)(cos(0) + l)LegendrcP(Z, 2, cos(0)) = (2) 
F 2 = HeunC(-2 icu, 2 icj, 4,-2 uj 2 , 4-1~1 2 + 2uj 2 , 1- \ r | ^((^OAt^M)^ =Qj 

where a; is a complex frequency, I is the angular momentum of the perturbation, 
9 £ [0, 7r] is the angle which we set to 9 = ir — 10~ 7 and \r\ — 20. HeunC is the 
confluent Heun function ( 3J) in maple notations. The parameter e is a small 
variation (|e|<l) in the phase condition arg(r)+arg(w) = — n/2 (see [T5]). 

Using Eqs. @, we run the two-dimensional Miiller algorithm to find the 
unknown I and u) with precision of the algorithms set to 15 digits. 

From the theory, it is known that I is an integer and I = 2, 3.... Comparing 
with the results obtained by the two-dimensional Miiller algorithm, for the first 
root I = 2, one has / = 1.99(9) + 1 x 10~ 17 i, with the first different from 9 digit 
being the 17 th . This shows that the new algorithm is capable of solving systems 
with one purely integer root in the pair with the expected precision. 

A comparison of the new algorithm with the well-known Newton's and Broy- 
den's methods, can be found in Table [T] 



Mode: 





1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


*fl W 


100 


99 


156 


196 


386 


240 


253 


282 


302 


368 


398 


tM2 [s] 


317 


413 


595 


741 


1175 


799 


874 


892 


1364 


971 


1355 


*M1 [s] 


202 


218 


335 


357 


497 


457 


396 


613 


623 


594 


667 



Table 1: The times needed for Broyden's method (ts) and the two-dimensional Miiller meth- 
ods (tjvfl and tMl) to fix a root. Note that while the precision of the former is 10 digits, 
the precision of the other two is 14 — 15 digits. To obtain those times, we solve the system: 
[Fx + F 2 , Fx - F 2 ] with starting points: u[n] + 0.01 + O.Oli, 2.1 + O.Oli, where n = 0..10. 
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Because the phase condition r =| r | e - I ( 1 + e ) 7r / 2 - M '"9(") includes the complex 
argument in non-analytical way, which cannot be differentiated, this problem 
cannot be solved directly using Newton's method. Broyden's method works 
but with serious limitation of its precision. This happens, because one of the 
roots y in the pair (x,y) is a real integer, while the other is complex and the 
algorithm fixes the integer root very quickly, thus the finite differences in the 
Jacobian become infinity. Because of this, the algorithm is able to fix only the 
first 10 — 11 digits, while the other algorithms fix 14 — 15 digits. Therefore, 
although Broyden's algorithm gives better times (see Table [I]) than the two- 
dimensional Muller algorithms, its precision is much lower and for modes with 
big imaginary part, it cannot be increased even by raising the software floating- 
point number to very high values. Furthermore, from Table [I] one can see that 
the time needed for the each algorithm to exit with a root dramatically increases 
with n. This emphasize on the importance of the convergence of the algorithm, 
which may become critical in physical problems where multiple roots must be 
found (see the second example). 

The numerical results for the QNMs are summed in Table In it, the 
QNM frequencies obtained from Sys. ^ are compared to those found by An- 
dersson ([E]) with the phase amplitude method. Recently, those results were 
confirmed by Fiziev (see 18J) with the one-dimensional Muller method applied 
on the exact solutions of the radial equation in terms of the confluent Heun 
function for 1 = 2. To check the accurateness of the new method, we evaluate 

A =| U!Muller2d ~ & Andersson \- 



n 


Our lu 


Andcrsson's u> 


A 





0.7473433689+0.177924631i* 


0.747343368+0. 177924630i 


1.68 x 10" 


9 


1 


0.6934219938+0.547829750i* 


0.693421994+0. 547829714i 


3.60 x 10" 


8 


2 


0.6021069092+0.956553966i* 


0.602106910+0. 956553966i 


1.02 x 10" 


9 


3 


0.5030099245+1.410296405i* 


0.503009924+1. 410296404i 


1.01 x 10" 


9 


4 


0.4150291596+1.893689782i* 


0.415029160+1. 893689782i 


4.41 x 10" 


10 


5 


0.3385988064+2.391216108i 


0.338598806+2. 391216108i 


9.67 x 10" 


10 


6 


0.2665046810+2.895821253i 


0.266504680+2. 895821252i 


1.48 x 10" 


9 


7 


0. 1856446684+3.407682345i 


0.185644672+3. 407682344i 


3.90 x 10" 


9 


8 


0.030649006+3.996823690i 


0+3.998000i** 


0.0306 




9 


0. 1265270180+4.605289542i 


0.126527010+4.605289530i 


1.44 x 10" 


8 


10 


0.1531069502+5. 121653272i 


0. 153106926+5. 121653234i 


4.52 x 10" 
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Table 2: A list of the frequencies we obtained for the QNMs of Schwarzschild black hole 
compared with the numbers found by Andersson. A =| i^Muller2d ~ u 'Andersson I ■ The first 5 
frequencies (n = — 4, marked with *) were obtained also by Fiziev using the confluent Heun 
functions and coincide with the presented here. The 8th mode, marked with **, was obtained 
by Leaver [25| . 



From the table, it is clear that in most cases, the modes obtained with the 
two-dimensional Muller algorithm coincide with those obtained by Andersson 
with more than 8 digits of precision in most cases and for modes n = 4, 5, 6, 
there are 9 coinciding digits (Andersson published 9 digits of his frequencies). 
These results also confirm the roots for n = 0, 1, 2, 3, 4 published in [18]. 

The mode with biggest deviation from the expected value is n = 8 in table [2] 
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and it was already discussed in [9] (and references therein) . In brief its properties 
are due to the branch cuts in the radial function, which also lead to non-trivial 
dependence of the frequencies on e (where arg(cj) + arg(r) = — 7r/2 and e < 1): 
for n < 4 u n = ±|SR(w„)| + for n > 4 w„(e) = -sgn(e)|5ft(u; n )| + S(w„)i 

and for n = 8, |e| < 0.75, w„ =8 = sgn(e) 0.030649006 + 3.996823690i. 

Because of this, the value for n = 8 in the table [2] was obtained for positive e 
(e = 0.3), unlike the other modes with n > 5, which were obtained for e = —0.3. 

The frequencies presented here are stable with precision of 6 digits at the 
worst and usually around 9 digits with respect to a change of e in the corre- 
sponding intervals. 

To confirm that the so-observed dependence of w n on the parameter e is not 
a problem of the algorithm, but it is a feature of the numerical realization of 
the confluent Hcun function, we make complex 3d plots of the radial function 
in use for several values of e. The effect of this parameter on the branch cuts 



in the radial function was already discussed in [9] and [TO], so on Fig. A. 2 A. 3 



|A.4| in the Appendix we only illustrate the movement of the branch cuts under 
the change of e. 

The appearance of those branch cuts represents one more complication in 
the work with the confluent Heun functions, since not all of them are docu- 
mented. Using the parameter e, however, those branch cuts can be controlled 
and one can obtain valuable physical results. 



3.2. Second example: Rotating black holes 

A significantly more complicated system to solve can be found in the case of 
QNMs from rotating black holes. In this case, one uses the exact solutions of the 
Teukolsky radial and angular equations, describing the linearized electromag- 
netic perturbations of the Kerr metric, in terms of confluent Heun functions, 
as stated for the first time in full detail in [IS]. The two-dimensional Miiller 
algorithm was applied successfully in this case too and the complete results and 
their discussion can be found in [10] . Here, we present some details on the 

numerical procedures used in this case. 

From [23], for the values of the parameters: s=-l, M=l/2, |r| = 110, m=0, 
a=0.01, 8 = 7r/3, one obtains: 

Fi(x,y) = HeunC(-1.9996ix,2.0002ja; + 1.000Q0.0002ix'-1.0000i-1.9996a;(i+x),1.9995x 2 -?/ 
+0.5000+1. 9998ix -110.02e( 4 ' 7124l - iar ' 9 ( :c »+l. 0000). (ii0.00e( 4 - 7124 ™ r9 ^»)( 2 ' 0W ' 0002l *> =0 
HeunC'(0.04x, -1.00, 1.00, -0.04x, 0.50-1. 00y+0.02a;-0.0001x 2 , 0.25) 
~~ HeunC(0.04x, -1.00, 1.00, -0.04x, 0.50- 1.00y+0.02x-0.0001x 2 , 0.25) + 

HcunC'(-0.04a:, 1.00, -1.00, 0.04x, 0.50- 1.00?/-0.02x-0.0001a; 2 , 0.75) _ 
HeunC(-0.04x, 1.00, -1.00, 0.04a;, 0.50-?/ - 0.02x-0.0001x 2 , 0.75) ~ 

where HeunC is the derivative of the confluent Heun function ([3]) as defined 

in MAPLE. 

For brevity, here the radial equation Fi(x, y) was rounded to only 4 digits of 
significance. In our numerical experiments, we used the complete system with 
software floating-point number set to 64, where the derivatives of the confluent 
Heun functions HeunC' were replaced with the associate Sn confluent Heun 
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function according to equation (3.7) of [7j. This was done to avoid the numerical 
evaluation HeunC' so that the peculiarities of the numerical implementation 
of the confluent Heun function (i.e. the use of maple fdiff procedure) are 
minimized. The difference in the times needed to fix a root when HeunC' is used 
and when it is not used is small for the modes (i.e.x) with small imaginary part 
(At ~ 15s), but it increases with the mode number, until it becomes significant 
for modes with big imaginary part (for the 10 th mode - R = 3 in table [3]- the 
difference is already At ~ 100s). This slowdown is due to the time-consuming 
numerical integration and differentiation in the complex domain, needed for the 
evaluation of HeunC'. 



R 




final 




*M2 [s] 


*M1 [s] 






(y\ final 


Broyden 


%2 


N U i 




0.49 + 0.18i 


0.4965436315+0. 1849695292i 


208 


102 


92 


1 


2.001 + 0.1i 


1.9999915063-0. 7347653. 10" 5 i 


23 


9(5)* 


11(4)* 




0.17 + 0.97i 


0.3495869222 + 1. 0503235984i 


449 


229 


244 


2 


2.001 + 0.1i 


2.0000392386 - 0.2937407.10~ 4 i 


34 


12(5)* 


15(5)* 




0.07 + 5.147i 


0.0608496029+5.1191008697i 


868 


568 


489 


3 


2.001 + 0.051i 


2.0010479243-0. 2491318. 10" 4 i 


36 


11(5)* 


17(5)* 



Table 3: QNMs of Kerr BH for s = —1. R numbers the root, t and N label the time and the 
iterations needed for the algorithms to exit. * denotes the roots dependent on the order of 
the equations in Ml and M2. 

For that system, three pairs of starting points were used:(0.49+0.18z, 2.001 + 
O.li), (0.17 + 0.97i, 2.001 +0.1i), (0.069 + 5.146i, 2.001 +0.05H). The results can 
be found in table [H One sees that the two modifications of the two-dimensional 
Miiller algorithm Ml and M2 are much quicker than the Broyden algorithm 
(tmi ~ tM2 < ts)- Newton's method once again cannot be used. 

The supremacy of the Miiller algorithms is clear and it is not isolated - it is 
observed for other modes or values of the parameters (for example, for m = 1). 
To check the precision of the method, the first two modes were compared with 
the already published results of electromagnetic QNMs of a Kerr black hole (see 
[20] ) and were found to coincide with at least 9 digits of significance with them. 
We could not find a published value for the third mode. 

This example shows that in systems featuring the confluent Heun functions 
and their derivatives, the two-dimensional Miiller method is much better suited 
than the already known algorithms. 

4. Conclusion 

The confluent Heun function appear in many physical problems. Because 
of the peculiarities of those functions, however, the standard numerical root- 
finding algorithms do not work efficiently enough on them. Here, we presented 
the general idea of a method for solving a system of two complex-valued nonlin- 
ear transcendental equations with complex roots based on the one-dimensional 
Miiller method. This method avoids some of the problems accompanying the 
work with the Heun functions and it is aimed to provide adequate way to deal 
with systems featuring those functions. 
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In the current article, the numerical results from the application of the new 
algorithm to systems from the QNM physics are presented. They showed that 
in those cases, the new method indeed works better than the standard methods. 
Therefore, the new method can be readily applied to find the roots of the Regge- 
Wheeler equation [TS] , the Zerilli equation [25J , the Teukolsky radial and angular 
equations [23], all of which are solved analytically in terms of confluent Heun 
functions. Using this algorithm, we were able to solve directly the problem 
of quasi- normal modes of a Schwarzschild ([9]) and Kerr black hole ([10]) with 
higher precision than that of the Broyden method. The so found solutions agree 
to great extent with previous published numerical results thus confirming the 
usefulness of the method. 

The complete mathematical investigation of the proposed new method, and 
especially its theoretical order of convergence under proper conditions on the 
class of functions Fi, Fi is still an open problem. 

For other applications of the method see the recent references [2§J [3D] ■ 
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Appendix A. The epsilon-method 






(d) e = 0.05 



Figure A. 2: 3d plots of the function F2, the solution of the RWE, in the complex interval 
u) = 0.32 + 1.4i..0.5 + 2Ai for e = 0.05,0,-0.05,-0.1 (the colors encode the phase of the 
complex function i"2.). The wall characteristic of the branching of the multivalued function 
is moved by e either to the left (e < 0) or to the right (e > 0). Note that on Figs. |A.2|A.3| 
and|A.4[ r =| r | e »((3/27r4<)/2-ar 9 (^) _ 




Figure A. 3: 3d plots of the function F2 for e = — 0.05,e = and e = 0.05 in the same interval 
for the complex uj as in |A.2| scaled near the expected roots. The different e lead to different 
profiles of the roots 
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ft(rc) »(») 



(c) e = (d) e = 0.05 

Figure A.4: Plot of the level curves of the function F 2 = \R(u, e)\ for e = 0.05, 0, -0.05, -0.08 
for [j = 0.32 + 1.4i..0.55 + 2.4i. The movement of the branch cut due to the change of e can 
be clearly seen. Note that on these plots, Fi is the radial solution of the Teukolsky Radial 
Equation [9], but as seen from |A.2| the branch cuts for this choice of parameters coincide with 
those of the solution of the RWE 
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